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We describe the implementation of total angular momentum dependent pseudopotentials in a 
plane wave formulation of density functional theory. Our approach thus goes beyond the scalar- 
relativistic approximation usually made in ab initio pseudopotential calculations and explicitly in- 
cludes spin-orbit coupling. We outline the necessary extensions and compare the results to available 
all-electron calculations and experimental data. 
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The ab initio pseudopotential method [ 0, has be- 
come a standard tool in many areas of electronic struc- 
ture calculation. Even magnetic compounds containing 
3d transition metal ions lie in the realm of the plane 
wave pseudopotential approach of density functional the- 
ory |q Hi- ln order to obtain high precision results it is 
necessary to include relativistic effects when calculating 
the electronic structure of materials containing third row 
elements || . Hence it is now standard procedure to cre- 
ate scalar-relativistic pseudopotentials that include the 
kinematic relativistic effects (mass-velocity and Darwin 
term) from the fully relativistic all-electron solution of 
the atom @J J [10. 

The spin-orbit interaction, however, is only effectively 
taken into account by the construction of j-averaged 
pseudopotentials for each angular momentum I. Thus no 
spin-orbit splitting is present in the resulting band struc- 
ture. Although the scalar-relativistic approximation is 
acceptable in many situations it becomes insufficient in 
cases where the observed quantities, such as hole effective 
masses or spin relaxation times, are a direct consequence 
of the spin-orbit splitting [Q . 

In this paper we report on the implementation of spin- 
orbit coupling in the pseudopotential scheme. We give 
the equations required to program the formalism within 
a generalized spinor approach, and compare the results to 
experimental data and to all-electron calculations. This 
is, to our knowledge, the first zeroth order implemen- 
tation of spin-orbit coupling in the ab initio pseudopo- 
tential scheme. Prior publications on this matter always 
relied on a second variation of the scalar-relativistic ze- 
roth order eigenstates, including spin-orbit coupling to 
first order in perturbation theory |]l3[ [l4[ |l5| . 

Although the fully relativistic treatment of the prob- 
lem would require a four-current formulation with Dirac 
spinors it has been shown by Kleinman that a Pauli- 
like Schrodinger equation captures all relativistic effects 
to order a 2 , where a is the fine structure constant |l6| |. 
The total ionic pseudopotential to be used is 
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where the are the total angular momentum eigen- 

functions which can be written in terms of the spher- 
ical harmonics, Y™, and the eigenfunctions of the z— 
component of the Pauli spin operator, | \) and | {). For 
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and for j = I — mj — m — \ have the form 
' I - m + 1 \ 3 



21 + 1 



Hence the operator Vps acts in both orbital and spin 
space. Note that there is only one radial pseudopotential 
component Vi.j with j = h for I = but two with j = l+h 
and j = I— j for each I > 0. The index m 3 in equation (fil) 
runs from — j to +j . It is computationally more efficient 
to transcribe each term of the semi-local pseudopotential 
operator Vps into the fully separable Kleinman-Bylander 
(KB) form 
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using the solutions of the atomistic pseudopotential prob- 
lem 
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where \R\ S j %a ) is the radial part of the pseudo eigenfunc- 
tion of atom species i s at position Ti s ^ a . The potential 
SV^y 1 " is defined as the difference 



8V{j ,la (r) = Vi d (r - r istia ) - Vi oc (r - r isAa ) 



(0) 



where V/ oc (r) is an arbitrary local potential that needs to 
be chosen such that the remaining SV's are short ranged. 
The complete KB pseudopotential operator is thus given 
as the sum of the local part and the non-local KB oper- 
ator. 

To our knowledge all previous pseudopotential calcu- 
lations that included spin-orbit coupling did so by using 
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a second variation step on the scalar-relativistic zeroth 
order wave functions, thus including the spin-orbit term 
to first order perturbation theory |l3|, |l5). In con- 
trast we solve directly for general two-component spinor 
Bloch wave functions expanding in a plane wave spinor 
basis 

|^ fc >=E + (7) 
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where G are reciprocal lattice vectors and er sums over 
up and down spin. In the basis of equation (Q) the action 
of the KB operator is as follows: 
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is a phase factor associated with the atomic position. 
The spin dependent factor M KB of equation (|s|) can be 
written as a spinor and for j = I + i is 
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for j = I — h. Also the last factor of equation (||), 
f?J a ,k,n,l,j, mj > depends on j as follows 
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for j = I — i Finally the KB factors appearing in 

equations (11) and (|lj) are defined as 
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S-FKKR a S-FLAPW 6 


R-PWPP 


Exp. c 


ao (A) 


5.56 5.620 


5.642 


5.653 


Bo (GPa) 


77 74 


72.2 


74.8 



TABLE I: Equilibrium lattice constant and bulk modulus de- 
termined in this work (R-PWPP) compared to all-electron 
calculations a Scalar relativistic FKKR, M. Asato et al., PRB 
60, 5202 (1999), 6 Scalar relativistic FLAPW, C. Filippi et al., 
PRB 50, 14947 (1994) and experiment c Landolt-Bornstein, 
Vol 22 (1987) 



where Y™ are the spherical harmonics, the polar angles 
9 and ip are determined by the vector k + G and ji are the 
spherical Bessel functions. The KB factors are calculated 
once and stored in memory. The contribution of a state 
\i>nk) to the non-local KB part of the total energy is thus 
given by 
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In order to test our implementation we calculated the 
properties of GaAs and compared with all-electron calcu- 
lations and available experimental data. Self-consistency 
was achieved by direct minimization of the total en- 
ergy via a conjugate gradient method The gallium 
and arsenic pseudopotentials were created following the 
Troullier-Martins scheme jl(| , and both contained s and 
p components. Care must be given to the local part of the 
pseudopotential entering in equation (^]) to ensure good 
transferability. We used the j-average of the unbound 4d 
state in case of gallium and likewise the j-average of the 
p states for the local part of the arsenic pseudopotential. 

The results for lattice constant and bulk modulus for 
GaAs are shown in tabic Q. The good agreement be- 
tween our fully relativistic pseudopotential results and 
the scalar-relativistic all-electron values in table [j] con- 
firms our approach and reaffirms the notion that spin- 
orbit splittings have little effect on the structural prop- 
erties of semiconductors |Q . Our calculated bandstruc- 
ture in figure [l] on the other hand shows clear evidence of 
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spin-orbit coupling. The top of the valence band splits 
into the light hole, heavy hole manifold and the split off 
band, separated by 350meV. A similar split is also ob- 
served in the upper conduction bands at the Brillouin 
zone center. In table [n] we compare the characteristic 
spin-orbit splittings we obtained for GaAs at the exper- 
imental lattice constant with values from two fully rela- 
tivistic all-electron calculations found in the literature. 
The agreement with both all-electron calculations is ex- 
cellent. 



Splitting 


R-FKKR a R-FLAPW 6 


R-PWPP 


Ao(r? 5 ) 


0.35 


0.34 


0.35 




0.20 




0.20 


Ai(L?) 


0.09 


0.09 


0.09 


A(X|) 


0.21 


0.20 


0.22 



TABLE II: Spin-orbit splittings for GaAs obtained in this 
work (R-PWPP) compared to the results of two relativistic 
all-electron calculations: a Fully relativistic FKKR, S. Bei der 
Kellen, A. J. Freeman, PRB 54, 11187 (1996) and 6 Scalar 
relativistic FLAPW + 2nd variation, C. Filippi, et al, PRB 
50, 14947 (1994) 



For completeness we compare in table III the eigen- 
value spectrum at three special k points of the Brillouin 
zone with the eigenvalues obtained from the two fully rel- 
ativistic all-electron calculations cited in table ||. The 
zero of energy was chosen to coincide with the top of 
the valence band. Despite the generally good agreement 
there are two obvious discrepancies at the Brillouin zone 
center that need some clarification. First the direct band 
gap of our pseudopotential calculation at T is more than 
4 times larger than the gap resulting from the all-electron 
calculations. Second the valence band width of our ap- 
proach is slightly smaller compared to the all-electron 
results. Both discrepancies result from the fact that the 
gallium 3d orbitals were placed in the frozen core in our 




Level 


R-FKKR" 


R-FLAPW 6 


R-PWPP 


1 6 


-12.94 


-12.91 


-12.67 


1 7 


-0.35 


-0.34 


-0.35 


1 8 


n nn 
u.uu 


n nn 
u.uu 


n nn 
u.uu 




0.12 


0.17 


0.69 


r? 


3.46 




3.40 


rg 


3.66 




3.60 


xi 


-10.42 


-10.41 


-10.35 


X? 


-7.02 


-7.00 


-6.83 


YV 


9 88 


9 8^ 


9 7 A 


xy 


-2.79 


-2.76 


-2.65 


xs 


1.17 


1.23 


1.23 


A? 


1.39 




1.44 


LI 


-11.18 


-11.14 


-11.06 


T v 


-6.83 


-6.82 


-6.63 


T V 

L 6 


-1.38 


-1.36 


-1.32 


T v 


-1.17 


-1.16 


-1.10 


L% 


0.71 


0.73 


0.97 


L% 


4.38 




4.34 


T c 


4.46 




4.44 



FIG. 1: Fully relativistic bandstructure of GaAs obtained 
with the implementation described in this work. 



TABLE III: Eigenvalue spectrum at three special k points 
compared with the same relativistic all-electron calculations 
cited in table [n| 



calculation but are free to change in the all-electron ap- 
proaches. Due to the well known self interaction problem 
of the local density approximation to density functional 
theory jl9| these fairly localized states will lie too high 
in energy when not frozen. The symmetry of the d states 
in the zincblende lattice at T only allows hybridization 
with p states, e.g. the top of the valence band. Hence 
the top of the valence band will shift upwards, leading to 
a reduced band gap and at the same time an increase in 
the valence band width. Due to the mixed character of 
the band states away from the Brillouin zone center the 
effect of the gallium 3d states is most pronounced at T. 

Compared to calculations that do not include the spin- 
orbit term we find that the inclusion of spin-orbit cou- 
pling worsens the short coming of the local density ap- 
proximation of underestimating the band gap. The rea- 
son for this observation simply lies in the fact that the 
top of the valence band splits and the light and heavy 
hole states move closer to the bottom of the conduction 
band. 

In conclusion, we have implemented spin-orbit cou- 
pling in the well established ab initio pseudopotential ap- 
proach of density functional theory. This paper gives the 
necessary expressions in a two-component spinor plane 
wave basis and demonstrates the applicability of the 
method for bulk GaAs. Our results compare very well 
to relativistic all-electron calculations. Since our direct 
approach is based on a complete spinor plane wave basis 
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it can easily be extended to systems that show exchange 
splitting and exhibit non-collinear spin arrangements. 

The code will be avail able under the GNU Gen- 
eral Public License |2Q| at ittp : //www.mrl .ucsb . edu/ 
"tin mr ich/ Spinor/ . 
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